clc;
clear;
close all;
%% 1. Inflation summary statistics and VAR
data  = xlsread('Data.xlsx',1);
% mean, standard deviation and autocorrelation
Y = data(:,2:5); % four inflation series
clear m sd ar1
[T,n] = size(Y);
for j = 1:n
m(j,:) = nanmean(Y(:,j)); % annualized
sd(j,:) = nanstd(Y(:,j)); % annualized
fit =   fitlm(Y(1:end-1,j), Y(2:end,j));
ar1(j,:) =fit.Coefficients{2,1:2}; % point estimate, 
ShocksAR(j,:) = fit.Residuals.Raw; 
end
inflation_stat = [m sd ar1];
inflation_corr = corr(Y);

% Share from regressions
fit = fitlm(Y(:,[2 3 4]),Y(:,1));
coef = zeros (4,2);
for i = 1:4
    coef(i,1) = fit.Coefficients.Estimate(i,1);
    coef(i,2) = fit.Coefficients.SE(i,1);
end
AA = [inflation_stat coef inflation_corr];
Table_1A = AA(:,1:3);
Table_1B = AA(2:end,5:6);
Table_1C = AA(:,7:10);

% inflation vol (Table 13)
inf_vol = zeros(2,4);
inf_vol(1,:) = std(Y(1:(1999-1964+1)*4+3,:)); % the first half, 1963Q2 - 1999Q4
inf_vol(2,:) = std(Y((1999-1964+1)*4+3+1:end,:)); % the second half, 2000Q1 - 2019Q4

% VAR
Y = data(2:end,2:8);
X = [data(1:end-1,2:8) ones(T-1,1)];
Phi = (X'*X)\(X'*Y);
U = Y-X*Phi;
ShocksVAR = U(:,1:4);
EXPINF = X*Phi;
 % ShocksVAR has four columns: CPI, core, food, energy
Shocks1 = ShocksVAR(:,1); % headline shock
Shocks2 = ShocksVAR(:,[2 4]); % core and energy shock
% Sticky/flexible inflation and comparison with core/energy inflation
price_sticky = data(:,9:10);
T_sticky = T-3*4-3;
YY2 = [data(4*3+3+1:end,[2 6 7 8 9 10])];
Y2 = YY2(2:end,:);
X2 = [YY2(1:end-1,:) ones(T_sticky-1,1)];
Phi2 = (X2'*X2)\X2'*Y2;
U2 = Y2-X2*Phi2;
Shocks_sticky = U2(:,[5 6]);
sticky_corr = zeros(4,1);
sticky_corr(1,1) = corr(Shocks_sticky(:,1),Shocks2((1966-1964+1)*4+2+1+1:end,2)); % energy and flexible
sticky_corr(2,1) = corr(Shocks_sticky(:,2),Shocks2((1966-1964+1)*4+2+1+1:end,1)); % core and sticky
sticky_corr(3,1) = corr(Shocks_sticky(:,2),Shocks2((1966-1964+1)*4+2+1+1:end,2)); % energy and sticky
sticky_corr(4,1) = corr(Shocks_sticky(:,1),Shocks2((1966-1964+1)*4+2+1+1:end,1)); % core and flexible
%% 2. Asset pricing test
Shock_1 = ShocksVAR(:,1);
Shock_2 = ShocksVAR(:,[2 4]);
asset_return_1 = xlsread('Data.xlsx',2);
asset_return_1 = asset_return_1(:,2:end);
asset_return_2 = xlsread('Data.xlsx',3);
asset_return_2 = asset_return_2(:,2:end);

[table_2a,~] = ap_test(Shock_1,asset_return_1,1);
[table_2b,~] = ap_test(Shock_2,asset_return_1,1);

[table_3a,~] = ap_test(Shock_1,asset_return_2,1);
[table_3b,~] = ap_test(Shock_2,asset_return_2,1);

[~,table_4a] = ap_test(Shock_1,asset_return_1,1);
[~,table_4b] = ap_test(Shock_2,asset_return_1,1);
[~,table_4c] = ap_test(Shock_1,asset_return_2,1);
[~,table_4d] = ap_test(Shock_2,asset_return_2,1);


%% Factor mimicking portfolio
within{1} = xlsread('Data',4);
within{2} = xlsread('Data',5);
within{3} = xlsread('Data',6);
within{4} = xlsread('Data',7);
within{5} = xlsread('Data',8);
within{6} = xlsread('Data',9);
within{7} = xlsread('Data',10);
within{8} = xlsread('Data',11);
within{9} = xlsread('Data',2);
within{10} = xlsread('Data',3);
table_5 = zeros(9,10);
T = length(Shock_1);
xx = ones(T,1);
for i = 1:10
    [~,R1] = fmp(Shock_1,within{i}(:,2:end));
    [~,R2] = fmp(Shock_2,within{i}(:,2:end));
    table_5(1,i) = nanmean(R2(:,1));
    EstCov1 = hac(xx,R2(:,1),'Intercept',false);
    table_5(2,i) = nanmean(R2(:,1))./sqrt(EstCov1);
    table_5(3,i) = nanmean(R2(:,1))./nanstd(R2(:,1))*sqrt(4);
    table_5(4,i) = nanmean(R2(:,2));
    EstCov2 = hac(xx,R2(:,2),'Intercept',false);
    table_5(5,i) = nanmean(R2(:,2))./sqrt(EstCov2);
    table_5(6,i) = nanmean(R2(:,2))./nanstd(R2(:,2))*sqrt(4);
    table_5(7,i) = nanmean(R1);
    EstCov3 = hac(xx,R1,'Intercept',false);
    table_5(8,i) = nanmean(R1)./sqrt(EstCov3);
    table_5(9,i) = nanmean(R1)./nanstd(R1)*sqrt(4);
end

%% Subsample analysis
[table_6a,~] = ap_test(Shock_1(1:(1999-1964+1)*4+2),asset_return_1(1:(1999-1964+1)*4+2,:),2);
[table_6b,~] = ap_test(Shock_2(1:(1999-1964+1)*4+2,:),asset_return_1(1:(1999-1964+1)*4+2,:),2);
[table_6c,~] = ap_test(Shock_1((1999-1964+1)*4+2+1:end),asset_return_1((1999-1964+1)*4+2+1:end,:),2);
[table_6d,~] = ap_test(Shock_2((1999-1964+1)*4+2+1:end,:),asset_return_1((1999-1964+1)*4+2+1:end,:),2);

[~,table_6e] = ap_test(Shock_1(1:(1999-1964+1)*4+2),asset_return_2(1:(1999-1964+1)*4+2,:),2);
[~,table_6f] = ap_test(Shock_2(1:(1999-1964+1)*4+2,:),asset_return_2(1:(1999-1964+1)*4+2,:),2);
[~,table_6g] = ap_test(Shock_1((1999-1964+1)*4+2+1:end),asset_return_2((1999-1964+1)*4+2+1:end,:),2);
[~,table_6h] = ap_test(Shock_2((1999-1964+1)*4+2+1:end,:),asset_return_2((1999-1964+1)*4+2+1:end,:),2);

% Test structural break
[~,N_asset] = size(asset_return_1);
dummy = zeros(T,1);
dummy((1999-1964+1)*4+2+1:end,1) = 1;
b_headline = zeros(N_asset,1);
b_core = zeros(N_asset,2);
for j = 1:N_asset
    result_1 = fitlm([Shock_1 dummy Shock_1.*dummy],asset_return_1(:,j));
    b_headline(j,1) = result_1.Coefficients.pValue(4);       
    result_2 = fitlm([Shock_2 dummy Shock_2.*dummy],asset_return_1(:,j));
    b_core(j,:) = result_2.Coefficients.pValue(5:6)';
end
table_6i = [b_headline b_core];

% Inflation-output correlation
gdp_growth = xlsread('Data',13);
gdp_growth = gdp_growth(2:end,2);
inflation = Y(:,[1 2 4]);
table_6j = zeros(3,2);
for i = 1:3
    table_6j(i,1) = corr(gdp_growth(1:(1999-1964+1)*4+2,1),inflation(1:(1999-1964+1)*4+2,i));
    table_6j(i,2) = corr(gdp_growth((1999-1964+1)*4+3:end,1),inflation((1999-1964+1)*4+3:end,i));
end
% Stock-bond correlation
stock = asset_return_1(:,1);
bond = asset_return_2(:,8);
sb_corr_1 = corr(stock(1:(1999-1964+1)*4+2),bond(1:(1999-1964+1)*4+2)) % Stock-bond correlation in two subsamples
sb_corr_2 = corr(stock((1999-1964+1)*4+2+1:end),bond((1999-1964+1)*4+2+1:end))
result1_stock = fitlm(Shocks2(1:(1999-1964+1)*4+2,:),asset_return_1(1:(1999-1964+1)*4+2,1));
b_stock_1 = result1_stock.Coefficients.Estimate(2:3,1);
result2_stock = fitlm(Shocks2((1999-1964+1)*4+2+1:end,:),asset_return_1((1999-1964+1)*4+2+1:end,1));
b_stock_2 = result2_stock.Coefficients.Estimate(2:3,1);

result1_bond = fitlm(Shocks2(1:(1999-1964+1)*4+2,:),asset_return_2(1:(1999-1964+1)*4+2,8));
b_bond_1 = result1_bond.Coefficients.Estimate(2:3,1);
result2_bond = fitlm(Shocks2((1999-1964+1)*4+2+1:end,:),asset_return_2((1999-1964+1)*4+2+1:end,8));
b_bond_2 = result2_bond.Coefficients.Estimate(2:3,1);

stock_hat_1 = Shocks2(1:(1999-1964+1)*4+2,:)*b_stock_1;
bond_hat_1 = Shocks2(1:(1999-1964+1)*4+2,:)*b_bond_1;
stock_hat_2 = Shocks2((1999-1964+1)*4+2+1:end,:)*b_stock_2;
bond_hat_2 = Shocks2((1999-1964+1)*4+2+1:end,:)*b_bond_2;

var_stock_1 = std(stock(1:(1999-1964+1)*4+2));
var_stock_2 = std(stock((1999-1964+1)*4+2+1:end));
var_bond_1 = std(bond(1:(1999-1964+1)*4+2));
var_bond_2 = std(bond((1999-1964+1)*4+2+1:end));

cov_sb_1 = cov(stock_hat_1,bond_hat_1);
cov_sb_2 = cov(stock_hat_2,bond_hat_2);
sbfit_corr_1 = cov_sb_1(1,2)/(var_stock_1*var_bond_1) % Fitted stock-bond correlation in two subsamples
sbfit_corr_2 = cov_sb_2(1,2)/(var_stock_2*var_bond_2)
%% Macroeconomic factor
macro = xlsread('Data',12);
Shock_macro{1} = macro(:,2); % consumption growth
Shock_macro{2} = macro(:,2:3); % nondurable and durable cons growth
Shock_macro{3} = macro(:,4); % IP growth
Shock_macro{4} = macro(:,5); % Payroll growth
Shock_macro{5} = macro(:,6); % Unemployment growth
Shock_macro{6} = macro(:,7:8); % HHL
Shock_macro{7} = macro(:,9); % Unfiltered cons growth
Shock_macro{8} = macro(:,10); % Capital share growth
Shock_macro{9} = macro(:,11); % Intermediary factor
Shock_macro{10} = macro(:,12); % CAPM

table_7 = nan(9,10);
for i = 1:10
    Shock = [Shock_2 Shock_macro{i}];
    [~,NN] = size(Shock);
    [~,table_7([(1:2*NN)';9],i)] = ap_test(Shock,asset_return_2,2);
end
%% Inflation predict growth
% 12 month inflation
growth = xlsread('Data',13);
growth = growth(:,2:end); % 3 growth series
inflation = data(:,2:5); % 4 inflation series
gset =  growth*400;  
IG = size(gset,2); % NO of regressand
IX = 2;

table_8a = nan(3,3);
table_8b = nan(3,5);
table_8c = nan(3,3);
table_8d = nan(3,5);
for ig = 1:IG 
    for ix = 1:IX
        % pick YY0 from gset and XX0
        YY0 = gset(:,ig); 
        XX1 = inflation(:,1); 
        XX2 = inflation(:,[2 4]); 
        for ik = [1 4]
        Y0_1 = nan(length(YY0),1);
        for t = ik:length(YY0)
            Y0(t,1) = sum(YY0(t-ik+1:t,1)); %t:t+ik
        end
        Y0 = Y0(1+ik :end,1);
        X1 =  XX1(1:end-ik,:);
        X2 =  XX2(1:end-ik,:)   ;
        fit1 = fitlm(X1,Y0);
        fit2 = fitlm(X2,Y0);

        [EstCov1,se1,coeff1] =hac( X1, Y0, 'display','off');
        [EstCov2,se2,coeff2] =hac( X2, Y0, 'display','off');

        t1 = coeff1./se1;
        t2 = coeff2./se2;

        if ik == 1
        table_8a(ig,:) = [coeff1(2) t1(2) fit1.Rsquared.Ordinary];  
        table_8b(ig,:) = [coeff2(2) t2(2) coeff2(3) t2(3) fit2.Rsquared.Ordinary];  
        end
        if ik == 4
        table_8c(ig,:) = [coeff1(2) t1(2) fit1.Rsquared.Ordinary];  
        table_8d(ig,:) = [coeff2(2) t2(2) coeff2(3) t2(3) fit2.Rsquared.Ordinary];         
        end
        end
    end
end  

%% Sticky and flexible inflation
Shocks_sticky = Shocks_sticky(:,[2 1]); % sticky, flexible
[table_10a,~] = ap_test(Shocks_sticky,asset_return_1((1966-1964+1)*4+2+1+1:end,:),2);
[~,table_10b] = ap_test(Shocks_sticky,asset_return_1((1966-1964+1)*4+2+1+1:end,:),2);
[~,table_10c] = ap_test(Shocks_sticky,asset_return_2((1966-1964+1)*4+2+1+1:end,:),2);
%% Time-varying lambda
spread = xlsread('Data',14);
spread = spread(:,2);
X = [Y spread];
[beta,lambda0,lambda1,tstatbeta,selambda0,selambda1,r2]=adrianetal(X,asset_return_2,[2 4],[8]);
lambda = [[lambda0;lambda1] , [lambda0;lambda1]./[selambda0;selambda1]];
% The first and third row correspond to the numbers mentioned in section 2.4.3.
%% TIPS
tips = xlsread('Data',15);
tips = tips(:,2);
fitlm(Shocks1,tips)
fitlm(Shocks2,tips)

%% announcement 
load data_announcement 
% data_regression (core shock,  headline shock, ffr rate change, negative
% two-year bond return, S&P return)

Xs = data_regression(:,[1 2]);
Ys = data_regression(:,[ 3 4 5]);
aaa = [];
for jj = 1:size(Xs,2)
    aa  = [];
    for j = 1:size(Ys,2)
        fit = fitlm(Xs(:,jj),Ys(:,j));
        aa(j,:) = [fit.Coefficients{2,[1 3]}];
    end    
    aaa = [aaa aa];
end

table9a = [ aaa(1,1:2) nan(1,2)
           nan(1,2) aaa(1,3:4) ] ;
table9b = [ aaa(2,1:2) nan(1,2)
           nan(1,2) aaa(2,3:4) ] ;
table9c = [ aaa(3,1:2) nan(1,6)
           nan(1,2) aaa(3,3:4) nan(1,4) ] ;

% bi multivariate   
    aaa  = [];
    for j = 1:size(Ys,2)
        fit = fitlm(Xs ,Ys(:,j));
        aaa(j,:) = [fit.Coefficients{2,[1 3]} fit.Coefficients{3,[1 3]}];
    end    

table9a(3,:)  = aaa(1,:);
table9b(3,:) = aaa(2,:); 
table9c(3,:)  = [  aaa(3,:) nan(1,4)] ;

% stock on ff
X = data_regression(:,[3 ]);
Y = data_regression(:, 5);
fit = fitlm(X,Y );
table9c(4,:)  = [ nan(1,4) fit.Coefficients{2,[1 3]} nan(1,2) ] ;

% stock on core and ff
X = data_regression(:,[1 3 ]);
Y = data_regression(:, 5);
fit = fitlm(X,Y );
table9c(5,:)  = [fit.Coefficients{2,[1 3]} nan(1,2) fit.Coefficients{3,[1 3]} nan(1,2) ] ;  

% stock on core and ff and ty
X = data_regression(:,[1 3 4]);
Y = data_regression(:, 5);
fit = fitlm(X,Y );
table9c(6,:)  = [fit.Coefficients{2,[1 3]} nan(1,2) fit.Coefficients{3,[1 3]} fit.Coefficients{4,[1 3]}  ] ;  
 
 